# -*- coding: utf-8 -*-
# @Time    : 2017/11/3 11:32
# @Author  : Deyu.Tian
# @Site    :
# @File    : config.py
# @Software: PyCharm Community Edition
from __future__ import print_function
import os
# import glob
#
# from mpl_toolkits.mplot3d import Axes3D#for matplotlib 2.0.x
import matplotlib.pyplot as plt#, matplotlib.cm as cm,  matplotlib.patches as mpatches
# from image2gif import writeGif
# from PIL import Image as PIL_Image
#
# import numpy as np

import gdal

import config
LandsatDir = config.LandsatDir

from util import *
from geotiff_tools import *

E = 1e-09#RuntimeWarning: divide by zero encountered in divide
from osgeo import osr, gdal
import numpy as np

def ascii2tif(fold, out):
    """
    ascii to geotiff convertion
    create the gdal output file as geotiff
    set the no data value
    set the geotransform
    numpy.genfromtxt('your file', numpy.int8) #looks like int from you example
    reshape your array to the shape you need
    write out the array.
    :param fname:
    :return:
    """
    txts = list_all_txts(fold)

    for txt in txts:
        # Set file vars
        output_file = "{}\\{}aaa.tif".format(out, txt[-11:-4])
        # Create gtif
        driver = gdal.GetDriverByName("GTiff")
        dst_ds = driver.Create(output_file, 321, 161, 1, gdal.GDT_CFloat32)
        raster = np.genfromtxt(txt, dtype=np.float32)
        raster = np.flipud(raster)
        print(raster.shape)

        # top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
        dst_ds.SetGeoTransform([60, 0.25, 0, 55, 0, 0.25])

        # set the reference info
        srs = osr.SpatialReference()
        srs.SetWellKnownGeogCS("WGS84")
        dst_ds.SetProjection(srs.ExportToWkt())

        # write the band
        dst_ds.GetRasterBand(1).WriteArray(raster)
        break

def valid_imshow_data(data):
    data = np.asarray(data)
    if data.ndim == 2:
        return True
    elif data.ndim == 3:
        if 3 <= data.shape[2] <= 4:
            return True
        else:
            print('The "data" has 3 dimensions but the last dimension '
                  'must have a length of 3 (RGB) or 4 (RGBA), not "{}".'
                  ''.format(data.shape[2]))
            return False
    else:
        print('To visualize an image the data must be 2 dimensional or '
              '3 dimensional, not "{}".'
              ''.format(data.ndim))
        return False


def ThresholdNDSI(imgId):
    """
    SET NDSI THRESHOLD
    :param NDSI_DATA:
    :return:
    """
    band_data = img2array("{}/{}/{}_SNOW_pan_clip.tif".format(config.LandsatDir, imgId, imgId))
    band_data = np.transpose(band_data, (1,2,0))
    band_data = stretch_N(band_data, 0, 100)
    band_data = band_data[:, :, ::-1]
    if not valid_imshow_data(band_data):
        exit(0)


    ndsi_data = img2array("{}/{}/{}_NDSI.tif".format(config.LandsatDir, imgId, imgId))
    ndsi_data[ndsi_data > 0.7] = 1


    plt.figure(figsize=(50, 50), dpi=300)

    ax1 = plt.subplot(131)
    ax1.set_title('Snow Bands')
    ax1.imshow(band_data[:, :, 1:4], cmap=plt.get_cmap('gist_ncar'))

    ax2 = plt.subplot(133)
    ax2.set_title('NDSI Data'.format(imgId))
    ax2.imshow(ndsi_data, cmap=plt.get_cmap('gray'))

    plt.show()
    #plt.savefig('{}/patchof240.png'.format(figDir))


def calcNDSI(imgId):
    """
    calculte NDSI by Green and SWIR1 band
    :param imgId:
    :return:
    """

    cmd = "gdal_calc.py -A {}/{}/{}_SNOW_pan_clip.tif --A_band=2 -B {}/{}/{}_SNOW_pan_clip.tif --B_band=4 " \
          "--type='Float32' --outfile={}/{}/{}_NDSI.tif " \
          "--calc='(A-B)/(A+B+{})'".format(LandsatDir, imgId, imgId, LandsatDir, imgId, imgId,
                                           LandsatDir, imgId, imgId, E)
    os.system(cmd)
    pass

def read_6SCC_NDSI(snowDir):
    metadata, ndsi = readNDSI("{}\\LC08_L1TP_150033_20160720_20170323_01_USER_NDSI.tif".format(snowDir))
    ndsi[ndsi >= 0.30] = 1
    ndsi[ndsi < 0.30] = 0
    print("max and min of ndsi:", np.max(ndsi), np.min(ndsi))

    # #SHOW
    plt.figure(figsize=(100, 100), dpi=100)
    ax2 = plt.subplot()
    ax2.set_title('NDSI 6S CC WINTER')
    ax2.imshow(ndsi, cmap=plt.get_cmap('gray'))
    plt.show()

    array2raster("{}\\FLAASH_summer_NDSI_Binary.tif".format(snowDir), metadata, ndsi)
    pass

def correctWinterNDSIwithSummers(snowDir):
    metadata1, summer_ndsi = readNDSI("{}\\FLAASH_summer_NDSI_Binary_30m.tif".format(snowDir))
    metadata2, winter_ndsi = readNDSI("{}\\6S_CC_winter_NDSI_Normed_Binary.tif".format(snowDir))
    print(summer_ndsi.shape, winter_ndsi.shape)
    winter_ndsi[summer_ndsi == 1] = 1
    # #SHOW
    plt.figure(figsize=(100, 100), dpi=100)
    ax2 = plt.subplot()
    ax2.set_title('NDSI 6S CC WINTER')
    ax2.imshow(winter_ndsi, cmap=plt.get_cmap('gray'))
    plt.show()

    array2raster("{}\\6S_CC_winter_NDSI_Normed_Binary_corrected.tif".format(snowDir), metadata2, winter_ndsi)
    pass


if __name__ == '__main__':
    #imgId = "LC08_L1TP_150033_20160227_20170329_01_T1"
    #imgId = "LC08_L1TP_150033_20160720_20170323_01_T1"
    #calcNDSI(imgId)

    #ThresholdNDSI(imgId)
    # ascfold = "P:\Aoi_paper\data\snow\Che\snowdepth-2016\\2016"
    # outfold = "P:\Aoi_paper\data\snow\Che\\tiffs"
    # ascii2tif(ascfold, outfold)

    correctWinterNDSIwithSummers(config.snowDir)